Dissipative Solitary Waves in Granular Crystals 
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We provide a quantitative characterization of dissipative effects in one- dimensional granular crys- 
tals. We use the propagation of highly nonlinear solitary waves as a diagnostic tool and develop 
optimization schemes that allow one to compute the relevant exponents and prefactors of the dis- 
sipative terms in the equations of motion. We thereby propose a quantitatively-accurate extension 
of the Hertzian model that encompasses dissipative effects via a discrete Laplacian of the veloci- 
ties. Experiments and computations with steel, brass, and polytetrafluoroethylene reveal a common 
dissipation exponent with a material-dependent prefactor. 
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Introduction. Since the advent of the famous Fermi- 
Pasta-Ulam model over fifty years ago, nonlinear oscil- 
lator chains have received a remarkable amount of at- 
tention in a wide range of physical settings [H. Areas 
of intense theoretical and experimental interest over the 
last decade include (but are not limited to) DNA double- 
strand dynamics in biophysics 0], coupled waveguide 
arrays in nonlinear optics [3], breathing oscillations in 
micromechanical cantilever arrays [H , and Bose- Einstein 
condensation in optical lattices in atomic physics [Bj]. 

Within this general theme of interplay between nonlin- 
earity and discreteness, one of the key subjects has been 
the study of one-dimensional (ID) granular materials, 
which consist of chains of interacting particles that start 
from point contact with each other and deform elasti- 
cally when compressed. In contrast to classically-studied 
disordered granular media, highly-packed granular lat- 
tices have negligible frictional and rotational dynamics, 
in favor of axial stress propagation [E 0]. The highly 
nonlinear dynamic response of such "crystals" has been 
the subject of considerable attention II, M H 0, 0, 0, 
H 0, 0, 0, EE MJOM- Additionally, 
granular crystals can be created from numerous material 
types and sizes, which makes their properties extremely 
tunable [E H EE] ■ This flexibility is valuable not only for 
basic studies of the underlying physics but also in poten- 
tial ap plic ations such as shock [24j and energy absorbing 
layers jl3l EE EE 0] , sound focusing devices (tunable 



acoustic lenses and delay lines), actuator s I25,l26|, s ound 



absorption layers, and sound scramblers [H, 12, 23[. 

While the standard Hertzian force model has been 
used extensively in most dynamical investigations and 
is now textbook material [8j, |27|, recent experimentally- 
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motivated investigations have illustrated the challenging 
need to include dissipation effects 24, 28, 2^, 3(1. Dissi- 
pative terms related to friction [3 if. plasticity [321 ] , visco- 
elasticity [23], and viscous drag [23^3(1 have been pro- 
posed to model particle collisions 6,0, 0]- However, 
none of these models captures both qualitatively and 
quantitatively the decay and wave shape of the highly 
nonlinear solitary waves observed experimentally. It is 
this important experimental and theoretical aspect of 
packed granular lattices that we aim to tackle in this Let- 
ter through the combination of modeling, numerical and 
physical experiments, and a detailed comparison thereof. 
Based on the earlier propositions of Rcfs. [E 0, 0, 0], 
we illustrate the prevalent nature of dissipation in the 
form of a discrete Laplacian in the velocities with uni- 
form exponent and a material-dependent prefactor. The 
broad interest of our findings results not only from their 
general nature for granular crystals of different materials 
but also from the significance of similar models in other 
fields, such as ID lattice turbulence 34 1. 

Experimental Setup. We assembled a monodisperse 
chain of N beads (here we report results for N = 70 
but we performed experiments for up to N = 188 with 
similar results) of different materials (see Table [XJ) with 
radius R = 2.38 mm in a horizontal setup (see Fig. QJi) 
composed of four-garolite rod stand. (To ensure con- 
tact between the particles, the guide was tilted at 4 de- 
grees.) To directly visualize the waves, we embedded cal- 
ibrated piezo sensors (RC ~ 10 3 / us, Piezo Systems Inc; 
see Fig. lb of Ref. 111!) inside selected particles, as de- 
scribed in Refs. 0, 0, . We generated solitary 
waves by impacting the chain with a striker (identical to 
the particles in the chain) launched along a ramp. We 
calculated the impact velocities Wi mp (in m/s) using a 
high-speed camera at the end of the ramp: vi,2 = 1.77, 
1)3,4 = 1.55, 115,6 = 1-40, «7,8 = 1.04, and 1)9,10 = 0.79. 

Model. We model a dissipative chain of N spherical 
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FIG. 1: (Color online) (a) Schematic diagram of the experi- 
mental setup, (b) Solitary wave decay in a chain composed 
of 70 steel particles impacted by a steel bead with «; mp = Vi . 
The (blue) solid curves correspond to the recordings for sen- 
sors placed in particles 9, 16, 24, 31, 40, 50, 56, and 63. 



Material 


m (g) 


E (GPa) 


V 


a 


7 


Steel 


0.45 


193 


0.30 


1.81 ±0.25 


-5.58 ±1.30 


PTFE 


0.123 


1.46 


0.46 


1.68 ±0.16 


-1.56 ±0.19 


Brass 


0.48 


103 


0.34 


1.85 ±0.13 


-6.84 ±0.66 



TABLE I: Material properties (mass m, elastic modulus E, 
and Poisson ratio v) for stainless steel [HI, HU, PTFE [ill . 
I37I . |38|], and brass [3^]. The last two columns present our 
best estimates, together with their standard deviation, of the 
dissipation coefficients (a, 7). 



beads as a ID lattice with Hertzian interactions 



fa = A +7* 



Jn+l 



(1) 



where s = sgn(<5„-(5„+i), A = EV2R/[3m (l - i/ 2 )], n e 
{1, . . . , N}, y n is the deviation of the nth bead from its 
equilibrium, S n = max{y n _i — y n , 0} for n £ {2, . . . , N}, 
5i = 0, Sn+i = max{j/Ar,0}, E is the Young's (elastic) 
modulus of the beads, v is their Poisson ratio, m is their 
mass, and R is their radius. The particle n = repre- 
sents the striker. Dissipation is incorporated by using 
a phcnomcnological force, of prefactor 7 < 0, between 
adjacent beads that depends on their relative velocities 
(in particular, on the "discrete Laplacian" in the veloci- 
ties), generalizing earlier models (with dissipation expo- 
nent a — 1 specifically) for dry granular matter [30| . 

In contrast to previous works that a priori assume that 
a = 1 (i.e., that model dissipation using a linear dashpot) 
[30| . we determine both a and 7 by directly comparing 
experimental and numerical results. The general coef- 
ficient a is thus a phcnomcnological parameter derived 
from the best fitting. We introduce the absolute value 
and the sign parameter s into ([T]) to ensure that genuine 
dissipation is guaranteed irrespective of the sign of the 
relative velocities between consecutive beads. The units 
of 7 would depend on the value of a and, accordingly, 
are more properly investigated in dynamic models that 
incorporate dissipation based on detailed measurements 
of restitutive losses that cannot currently be achieved ex- 
perimentally 0]. Importantly, the value we obtain for a 
differs decidedly from the coefficients used in previous 
modeling attempts [3(j (see the discussion below). 

Determining the dissipation coefficients. We now deter- 
mine the "optimal" dissipation coefficients (a, 7) from 
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FIG. 2: (Color online) Optimization of the dissipation co- 
efficients (a, 7) for a steel chain, (a) Difference D(a,7), as 
defined in Eq. ([2]), between the force maxima recorded in the 
experiment and our model, (b) Difference A„(a,7), as de- 
fined in Eq. l|3|l. in wave forms between the experiment and 
our model for sensor n = 56. The solid and dashed curves 
correspond to the minima obtained from panels (a) and (b), 
respectively, (c) Maximum force F m (n) for experiments with 
t'imp = V3 (top curves) and Ui mp = v& (bottom curves, dis- 
placed by 5 units for clarity). The (red) circles correspond 
to the experiment, and the (green) thick curves give the nu- 
merical best fit with (0,7) = (1.81 ±0.25, -5.58 ±1.30). The 
dashed curves correspond to the extreme cases using the stan- 
dard deviation found in the optimal parameters, (d) Veloc- 
ity of traveling front versus the maximum force (in a log-log 
plot). The solid curve represents the best linear fit, which 
gives v oc F^ 17 ; we also show a dashed line with slope 1/6. 



the experimental data for different materials and differ- 
ent configurations. The experimental data consists of the 
time series of the force through each sensor. We optimize 
the pair (a, 7) by minimizing the following two differences 
between numerics and each particular experiment: 
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An(a,7)=77; / ^7Z\ dt - ( 3 ) 



T 



F ex P(n) 



where F%*> = (1/N) En=i F™ p (n), F c ^(n) = 
J t f F oxp (t;n)dt, F(t;n) is the time series data of the 
force through the nth sensor (see Fig. [lb) , and F m (n) = 
m&x t {F(t; n)} is the maximum force recorded by the nth 
sensor over the recording time span [ti, t/ = i,-+T], where 
T is typically about 100 us. The superscripts 'exp' and 
'num' denote, respectively, the experimental and numer- 
ical data. The function D(a,j) measures the "distance" 
between the numerics and the experiment using the max- 
ima of the forces through all sensors of the experiment. 
The function A n (a,"f) measures the difference between 
experimental and numerical pulse shapes that go through 
the nth sensor. In order to avoid biasing A n (cv,7) with 
the difference in force magnitude [which is already taken 
into account when optimizing D(a,"/)] 7 we rcscale the 
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FIG. 3: (Color online) Force versus time for the steel chain 
with iiimp = V2 through sensors at positions (a) n — 16 and 
(b) n = 56. The (red) thick solid curve depicts the (smoothed; 
see text) experimental series, and the thin (blue) dashed, dot- 
ted, and solid curves respectively show the numerics with 
(a, 7) = (1, -5.5), (1.4, -6), and (1.81, -5.58). The last case 
corresponds to the best fit (see text) for the dissipation pa- 
rameters for the chains of steel beads. 



experimental data so that the numerical and experimen- 
tal maxima match before we compare wave forms. That 
is, F c *P{t;n) -> F 0X P(t;n) x F^ UTa {n)/ F^P{n). Panels 
(a) and (b) in Fig. [2] depict, respectively, the differences 
D(a, 7) and A. n (a, 7) in a particular (a, 7) range for a 
steel chain using a sensor placed towards the end of the 
chain. As can be observed from these panels, the opti- 
mization of the force maxima F m [panel (a)] and the force 
pulse shape [panel (b)] are not sufficient on their own 
to determine the dissipation parameters. However, it is 
meaningful (and always well-defined) to optimize force 
maxima and pulse shape together by taking the intersec- 
tion between the minima of each case (see the point at the 
intersection of the solid and dashed curves). For experi- 
ment j (with impact velocity vj), we average the parame- 
ter pair (aj,7j) over four sensors located throughout the 

bead chain. Finally, we average (a, 7) = Ylj'=i( a j>^j) 
over the N e = 10 different experiments to obtain the 
optimal dissipation parameters (a, 7) and compute the 
standard deviation for the N e experiments. 

We summarize our results, for three different set 
of experiments — using steel, teflon (polytetrafluorocthy- 
lene; PTFE), and brass beads — in the last two columns 
of Table |TJ In order to validate the results of the above 
optimization procedure a posteriori, we take the opti- 
mal dissipation parameters for the steel bead chain and 
compare the maximal forces obtained numerically with 
the experiments in panel (c) of Fig. [2] In the panel, we 
show two typical examples (for impact velocities V3 and 
Vs) and also plot the curves incorporating the standard 
deviation measured in our analysis. As can be clearly 
observed, all experimental data points fall well within 
the predicted region. To further validate our results, we 
compared the dependence of the pulse velocity v against 
the maximal force F m in panel (d) of Fig. [2] [this panels 
shows a typical example; we obtained similar results for 
the other configurations (results not shown here)]. The 
obtained exponent is extremely close to the theoretical 
value of 1/6 (shown by the dashed line) § . 

In order to gain a deeper understanding of the role of 
the dissipation exponent a, we depict in Fig. [3] the pulse 



shape for two sensors in the steel chain (one near the 
beginning of the chain and the other one near the end). 
We depict the experimental pulse (smoothed by nearest- 
neighbor averaging) with the (red) solid curve. The thin 
(blue) curves, depict three numerical runs using three 
different pairs (a, 7) along the minimum curve [shown by 
a solid curve in panel (a) of Fig. [2]. 

It is interesting to note that for all materials tested, 
higher impact velocities correspond to a faster initial de- 
cay as compared to the latter part of the chain (proba- 
bly related to the initiation of plasticity at the contact). 
Also, by comparing the wave decay in chains composed 
of steel and teflon (or brass) beads, a faster and more 
pronounced energy loss is evident for the softer beads. 
To understand physically this dissipation, one should ex- 
plore a more detailed analysis of the contact plasticity, 
inelastic restitution, and hydrodynamic drag. We stress 
here that our granular crystals have a closely packed par- 
ticle arrangement and limited or null rotational dynamics 
and particles' displacements (small frictional response). 

Note that the optimal dissipation exponent a for the 
three material types considered is consonant with a value 
close to a = 1.75. This indicates the prevalence of the 
phenomcnological damping introduced in Eq. {T]), which 
is one of the principal findings of this Letter. It is impor- 
tant to point out the disparity of this optimal exponent 
from earlier investigations, which focused on the (linear 
dashpot) case of a = 1 [24, 3^, 34 1. On the other hand, 
naturally, the dissipation prefactor 7 does depend on the 
material. For steel and brass, which have similar mate- 
rial properties, 7 is also similar (steel has 7 = —5.58 and 
brass has 7 = —6.84). However, for teflon, as can be 
anticipated from the much lower elastic modulus E, the 
prefactor 7 is significantly smaller (7 = —1.56). 

We show typical examples of the results for teflon (left 
column) and brass (right column) in Fig. |4] The top pan- 
els depict the maximal force through the chain using the 
optimal dissipation parameters. Note in the pulse shape 
results (bottom panels) for teflon and brass that low dis- 
sipation exponents a tend to overestimate the size of the 
secondary pulse hump. Another relevant observation, in 
connection with its much smaller dissipation prefactor 7, 
is that chains of teflon beads may offer the first unam- 
biguous observation of the secondary pulses (see Fig. 0£) 
argued to arise for weaker dissipation in Ref . [30l | . 

Conclusions. In this Letter, we have offered for the first 
time a quantitative and systematic modeling attempt at 
the role of dissipation in granular crystals. Through de- 
tailed comparison of numerical simulations and experi- 
ments in a variety of materials (steel, teflon, and brass), 
we have demonstrated a generic functional form of the 
dissipation, modeled by a phenomcnological term based 
on the second difference of the velocities between adja- 
cent beads (i.e., a discrete Laplacian) that is raised to 
a common exponent. This allowed us to augment the 
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FIG. 4: (Color online) Results for the teflon (left column) 
and brass (right column) experiments. The top panels depict 
the same information as panel (c) in Fig. [2] for experiments 
with impact velocities (a) v$ (top curves) and vs (bottom 
curves, displaced by 0.3 units for clarity); and (b) velocities 
V3 (top curves) and ve (bottom curves, displaced by 7 units for 
clarity). The best fit for the dissipation parameters for teflon 
and brass are, respectively, (a, 7) = (1.68±0.16, — 1.56±0.19) 
and (a, 7) = (1.85 ± 0.13, -6.84 ± 0.66). The bottom panels 
show the same information as in Fig. [3] In panel (c), we 
depict the force versus time through the sensor at n = 38 
with (a, 7) = (1,-1.56), (1.4,-1.56), and (1.68,-1.56). In 
panel (d), we show the same information for the sensor at 
n = 14 with (a, 7) = (1, -5.5), (1.4, -6), and (1.85, -6.84). 



standard dynamical model based on Hertzian forces to 
encompass this dissipation effect in (optimal) quantita- 
tive agreement with our experiments. We found that the 
dissipation prefactor is material-dependent and that the 
considerably weaker prefactor of teflon (in comparison 
to brass and steel) allows one to observe unambiguously 
(and for the first time) secondary pulses such as the ones 
proposed in Ref. [3(| • Our study also provides a starting 
point for a potential first-principles derivation, as well 
as for future quantitative investigations of this newly- 
proposed model. For example, it would be worth exam- 
ining the critical prefactor below which a secondary wave 
should be expected to emerge, the interplay of the role of 
dissipation and plasticity (and a quantitative incorpora- 
tion of the latter) in the dynamics, and extensions of the 
present considerations to higher-dimensional settings. 
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